import sys
sys.path = ['/home/as259691/PycharmProjects/FluidDyn1D'] + sys.path
from src.main import *
from src.plot_fields import *
# %matplotlib inline
rc('figure', figsize=(10,5))
rc('figure', dpi=100)
Ici on va réaliser une simulation sans diffusion pour différentes écritures de notre équation thermique.
La résolution se fait à chaque fois en WENO avec Euler explicite en temps.
# d = 6./100*Delta/2.
phy_prop = PhysicalProperties(Delta=0.02, v=0.2, dS=0.005**2,
lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=1.,
alpha=0.06, a_i=357.)
markers = Bulles(phy_prop=phy_prop)
Formulation = [Problem, ProblemConserv2]
n = 1000
t_fin = 0.2
def compare_energy_forme(formu, phy_prop, num_prop, markers, t_fin):
fig1,ax1 = plt.subplots(1)
ax1.set_title('Énergie en fonction du temps')
for form in formu:
print()
prob = form(get_T_creneau, markers=markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob.energy
print(prob.name)
print('==========================')
t, e = prob.timestep(t_fin=t_fin, number_of_plots=5, debug=False, plotter=Plotter('decale'))
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob.dt / E0 # on a mult
# par Dt / rho_cp_l T_l V
print('dE*/dt* = %g' % dedt_adim)
le = fig1.legend()
En fait s'il n'y a pas de convection il n'y a pas de différence entre les différentes formes, à l'exception de la moyenne utilisée pour $\frac{1}{\rho C_p}$
num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno', time_scheme='rk4', phy_prop=phy_prop, cfl=0.5)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 6.918433404737903e-06 Cas : mixte, rk4, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = -0.000006 dt fourier 6.918433404737903e-06 Forme conservative boniou, Cas : mixte, rk4, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = 0.000022
num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 6.918433404737903e-06 Cas : mixte, euler, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = -0.000005 dt fourier 6.918433404737903e-06 Forme conservative boniou, Cas : mixte, euler, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = 0.000023
num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 6.918433404737903e-06 Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = -0.000006 dt fourier 6.918433404737903e-06 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = 0.000005
dx_list = [3.9*10**-6, 3.9*10**-5, 1*10**-5]
t_fin = 0.01
for dx in dx_list:
num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 6.896863867106643e-08 Cas : mixte, euler, weno upwind, dx = 3.90016e-06, dt = 6.89686e-08, cfl = 0.00353671 ========================== dE*/dt* = -0.000000 dt fourier 6.896863867106643e-08 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90016e-06, dt = 6.89686e-08, cfl = 0.00353671 ========================== dE*/dt* = 0.000000 dt fourier 6.918433404737903e-06 Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = -0.000103 dt fourier 6.918433404737903e-06 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = 0.000034 dt fourier 4.538601983461999e-07 Cas : mixte, euler, weno upwind, dx = 1.0005e-05, dt = 4.5386e-07, cfl = 0.00907267 ========================== dE*/dt* = -0.000004 dt fourier 4.538601983461999e-07 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 1.0005e-05, dt = 4.5386e-07, cfl = 0.00907267 ========================== dE*/dt* = 0.000001
dx_list = [7*10**-5]
t_fin = 0.01
for dx in dx_list:
num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 2.232841866976439e-05 Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05, cfl = 0.063636 ========================== dE*/dt* = -0.000469 dt fourier 2.232841866976439e-05 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05, cfl = 0.063636 ========================== dE*/dt* = 0.000172
Ici on ne change pas $\Delta x$, mais on diminue dt_min pour qu'il soit contraignant
dt_list = [2*10**-5, 1*10**-5, 5*10**-6]
t_fin = 0.01
for dt in dt_list:
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop, dt=dt)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt min 2e-05 Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2e-05, cfl = 0.057 ========================== dE*/dt* = -0.000418 dt min 2e-05 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2e-05, cfl = 0.057 ========================== dE*/dt* = 0.000154 dt min 1e-05 Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 1e-05, cfl = 0.0285 ========================== dE*/dt* = -0.000209 dt min 1e-05 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 1e-05, cfl = 0.0285 ========================== dE*/dt* = 0.000076 dt min 4.9999999999999996e-06 Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 5e-06, cfl = 0.01425 ========================== dE*/dt* = -0.000104 dt min 4.9999999999999996e-06 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 5e-06, cfl = 0.01425 ========================== dE*/dt* = 0.000038
t_fin = 0.01
phy_prop = PhysicalProperties(Delta=0.02, v=0., dS=0.005**2,
lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=1.,
alpha=0.06, a_i=357.)
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 2.232841866976439e-05 Cas : diffusion, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05 ========================== dE*/dt* = 0.000010 dt fourier 2.232841866976439e-05 Forme conservative boniou, Cas : diffusion, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05 ========================== dE*/dt* = -0.000000
t_fin = 0.01
phy_prop = PhysicalProperties(Delta=0.02, v=0.2, dS=0.005**2,
lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=0.,
alpha=0.06, a_i=357.)
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 2.232841866976439e-05 Cas : convection, euler, weno upwind, dx = 7.01754e-05, cfl = 0.063636 ========================== dE*/dt* = -0.000576 dt fourier 2.232841866976439e-05 Forme conservative boniou, Cas : convection, euler, weno upwind, dx = 7.01754e-05, cfl = 0.063636 ========================== dE*/dt* = 0.000166
Schemas = ['upwind', 'center', 'weno', 'weno upwind']
Time_scheme = ['euler', 'rk4']
def compare_energy_schema(schemas, form, time_scheme, phy_prop, markers, t_fin):
fig1,ax1 = plt.subplots(1)
ax1.set_title('Énergie en fonction du temps')
for schem in schemas:
for ts in time_scheme:
a = Plotter('decale')
num_prop = NumericalProperties(dx=3*10**-5, schema=schem, time_scheme=ts, phy_prop=phy_prop)
print()
prob = form(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob.energy
print(prob.name)
print('==========================')
t, e = prob.timestep(t_fin=t_fin, number_of_plots=5, debug=False, plotter=a)
a.ax.set_xlim(0., phy_prop.Delta/2)
a.ax.set_ylim(0.6, 1.1)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob.dt / E0 # on a mult
# par Dt / rho_cp_l T_l V
print('dE*/dt* = %f' % dedt_adim)
le = fig1.legend()
t_fin = 0.01
fig1, ax1 = plt.subplots(1)
a = Plotter('decale')
num_prop = NumericalProperties(dx=1*10**-5, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
a.ax.set_ylim(0.5,1.05)
a.ax.set_xlim(0., phy_prop.Delta/4.)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
ax1.legend()
ax1.grid(b=True)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0 # on a mult
num_prop = NumericalProperties(dx=5*10**-6, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
a.ax.set_ylim(0.5,1.05)
a.ax.set_xlim(0., phy_prop.Delta/4.)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
ax1.legend()
ax1.grid(b=True)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0 # on a mult
# par Dt / rho_cp_l T_l V
print('dE*/dt* = %f' % dedt_adim)
dt fourier 4.538601983461999e-07 Cas : convection, rk4, weno, dx = 1.0005e-05, cfl = 0.00907267 ========================== dt fourier 1.1335161290322582e-07 Cas : convection, rk4, weno, dx = 5e-06, cfl = 0.00453406 ========================== dE*/dt* = -0.000000
t_fin = 0.01
Schemas = ['upwind', 'weno', 'weno upwind']
compare_energy_schema(Schemas, Problem, Time_scheme, phy_prop, markers, t_fin)
dt fourier 4.0888316389624144e-06 Cas : convection, euler, upwind, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = -0.000129 dt fourier 4.0888316389624144e-06 Cas : convection, rk4, upwind, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = -0.000130 dt fourier 4.0888316389624144e-06 Cas : convection, euler, weno, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = -0.000009 dt fourier 4.0888316389624144e-06 Cas : convection, rk4, weno, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = -0.000017 dt fourier 4.0888316389624144e-06 Cas : convection, euler, weno upwind, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = -0.000081 dt fourier 4.0888316389624144e-06 Cas : convection, rk4, weno upwind, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = -0.000083
Schemas = ['upwind', 'weno', 'weno upwind']
Time_scheme = ['euler', 'rk4']
compare_energy_schema(Schemas, ProblemConserv2, Time_scheme, phy_prop, markers, t_fin)
dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : convection, euler, upwind, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = 0.000033 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : convection, rk4, upwind, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = 0.000032 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : convection, euler, weno, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = 0.000011 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : convection, rk4, weno, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = 0.000015 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : convection, euler, weno upwind, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = 0.000014 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 3.003e-05, cfl = 0.0272316 ========================== dE*/dt* = 0.000014
from line_profiler import LineProfiler
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='rk4', phy_prop=phy_prop, dt=1.)
lp = LineProfiler()
func = [Problem.rk4_timestep, Bulles.indicatrice_liquide, interpolate_form_center_to_face_weno,
integrale_volume_div]
for fu in func:
lp.add_function(fu) # add additional function to profile
lp_wrapper = lp(compare_energy_forme)
lp_wrapper(Formulation, phy_prop, num_prop, markers, 0.1)
lp.print_stats()
dt fourier
2.232841866976439e-05
Cas : convection, rk4, weno upwind, dx = 7.01754e-05, cfl = 0.063636
==========================
dE*/dt* = -0.000018
dt fourier
2.232841866976439e-05
Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 7.01754e-05, cfl = 0.063636
==========================
dE*/dt* = 0.000018
Timer unit: 1e-06 s
Total time: 13.7396 s
File: /home/as259691/PycharmProjects/FluidDyn1D/src/main.py
Function: integrale_volume_div at line 10
Line # Hits Time Per Hit % Time Line Contents
==============================================================
10 def integrale_volume_div(center_value, face_value, I=None, cl=1, dS=1., schema='weno'):
11 """
12 Calcule le delta de convection aux bords des cellules
13 :param dS:
14 :param schema:
15 :param center_value: les valeurs présentes aux centres des cellules
16 :param face_value: les valeurs présentes aux faces des cellules
17 :param cl: si cl = 1, on prend des gradients nuls aux bords du domaine
18 :return: le delta au centre
19 """
20 89560 97023.0 1.1 0.7 if len(center_value.shape) != 1 or len(face_value.shape) != 1:
21 raise NotImplementedError
22 89560 64403.0 0.7 0.5 if center_value.shape[0] != face_value.shape[0] - 1:
23 print('Le tableau des valeurs aux faces ne correspond pas au tableau des valeurs au centre')
24 print(center_value.shape, face_value.shape)
25 raise NotImplementedError
26 89560 35334.0 0.4 0.3 if schema is 'upwind':
27 interpolated_value = interpolate_from_center_to_face_upwind(center_value, cl)
28 89560 31055.0 0.3 0.2 elif schema is 'center':
29 interpolated_value = interpolate_from_center_to_face_center(center_value, cl)
30 89560 29909.0 0.3 0.2 elif schema is 'weno':
31 interpolated_value = interpolate_form_center_to_face_weno(center_value, cl=1)
32 89560 40326.0 0.5 0.3 elif schema == 'weno upwind':
33 89560 30846.0 0.3 0.2 if I is None:
34 raise NotImplementedError
35 89560 13061058.0 145.8 95.1 interpolated_value = interpolate_center_value_weno_to_face_upwind_interface(center_value, I, cl=1)
36 else:
37 raise NotImplementedError
38 89560 106096.0 1.2 0.8 flux = interpolated_value * face_value
39 # plt.figure(1)
40 # plt.plot(flux)
41 # plt.show(block=False)
42 89560 137855.0 1.5 1.0 delta_center_value = flux[1:] - flux[:-1]
43 89560 105728.0 1.2 0.8 return dS * delta_center_value
Total time: 7.54102 s
File: /home/as259691/PycharmProjects/FluidDyn1D/src/main.py
Function: interpolate_form_center_to_face_weno at line 73
Line # Hits Time Per Hit % Time Line Contents
==============================================================
73 def interpolate_form_center_to_face_weno(a, cl=1):
74 """
75 Weno scheme
76 :param a: the scalar value at the center of the cell
77 :param cl: conditions aux limites, cl = 1: périodicité
78 :return: les valeurs interpolées aux faces de la face -1/2 à la face n+1/2
79 """
80 89560 113876.0 1.3 1.5 if cl == 1:
81 89560 229841.0 2.6 3.0 center_values = np.empty(a.size + 5)
82 89560 217261.0 2.4 2.9 center_values[:3] = a[-3:]
83 89560 147485.0 1.6 2.0 center_values[3:-2] = a
84 89560 153844.0 1.7 2.0 center_values[-2:] = a[:2]
85 else:
86 raise NotImplementedError
87 89560 120102.0 1.3 1.6 ujm2 = center_values[:-4]
88 89560 116337.0 1.3 1.5 ujm1 = center_values[1:-3]
89 89560 116051.0 1.3 1.5 uj = center_values[2:-2]
90 89560 118486.0 1.3 1.6 ujp1 = center_values[3:-1]
91 89560 117561.0 1.3 1.6 ujp2 = center_values[4:]
92 89560 473192.0 5.3 6.3 f1 = 1. / 3 * ujm2 - 7. / 6 * ujm1 + 11. / 6 * uj
93 89560 376808.0 4.2 5.0 f2 = -1. / 6 * ujm1 + 5. / 6 * uj + 1. / 3 * ujp1
94 89560 369830.0 4.1 4.9 f3 = 1. / 3 * uj + 5. / 6 * ujp1 - 1. / 6 * ujp2
95 89560 150119.0 1.7 2.0 eps = np.array(10. ** -6)
96 89560 987023.0 11.0 13.1 b1 = 13. / 12 * (ujm2 - 2 * ujm1 + uj) ** 2 + 1. / 4 * (ujm2 - 4 * ujm1 + 3 * uj) ** 2
97 89560 657626.0 7.3 8.7 b2 = 13. / 12 * (ujm1 - 2 * uj + ujp1) ** 2 + 1. / 4 * (ujm1 - ujp1) ** 2
98 89560 850713.0 9.5 11.3 b3 = 13. / 12 * (uj - 2 * ujp1 + ujp2) ** 2 + 1. / 4 * (3 * uj - 4 * ujp1 + ujp2) ** 2
99 89560 334858.0 3.7 4.4 w1 = 1. / 10 / (eps + b1) ** 2
100 89560 309056.0 3.5 4.1 w2 = 3. / 5 / (eps + b2) ** 2
101 89560 304737.0 3.4 4.0 w3 = 3. / 10 / (eps + b3) ** 2
102 89560 207644.0 2.3 2.8 sum_w = w1 + w2 + w3
103 89560 254552.0 2.8 3.4 w1 /= sum_w
104 89560 180639.0 2.0 2.4 w2 /= sum_w
105 89560 171257.0 1.9 2.3 w3 /= sum_w
106 89560 358700.0 4.0 4.8 interpolated_value = f1 * w1 + f2 * w2 + f3 * w3
107 89560 103418.0 1.2 1.4 return interpolated_value
Total time: 4.6227 s
File: /home/as259691/PycharmProjects/FluidDyn1D/src/main.py
Function: indicatrice_liquide at line 245
Line # Hits Time Per Hit % Time Line Contents
==============================================================
245 def indicatrice_liquide(self, x):
246 """
247 Calcule l'indicatrice qui correspond au liquide avec les marqueurs selon la grille x
248 :param x: les positions des centres des mailles
249 :return: l'indicatrice
250 """
251 44784 418997.0 9.4 9.1 i = np.ones_like(x)
252 44784 56580.0 1.3 1.2 dx = x[1] - x[0]
253 223920 220913.0 1.0 4.8 for markers in self.markers:
254 179136 131025.0 0.7 2.8 if markers[0] < markers[1]:
255 176447 778653.0 4.4 16.8 i[(x > markers[0]) & (x < markers[1])] = 0.
256 else:
257 2689 12116.0 4.5 0.3 i[(x > markers[0]) | (x < markers[1])] = 0.
258 179136 687659.0 3.8 14.9 diph0 = (np.abs(x - markers[0]) < dx / 2.)
259 179136 834643.0 4.7 18.1 i[diph0] = (markers[0] - x[diph0]) / dx + 1. / 2.
260 179136 630356.0 3.5 13.6 diph1 = (np.abs(x - markers[1]) < dx / 2.)
261 179136 833671.0 4.7 18.0 i[diph1] = -(markers[1] - x[diph1]) / dx + 1 / 2.
262 44784 18085.0 0.4 0.4 return i
Total time: 10.1102 s
File: /home/as259691/PycharmProjects/FluidDyn1D/src/main.py
Function: rk4_timestep at line 528
Line # Hits Time Per Hit % Time Line Contents
==============================================================
528 def rk4_timestep(self, debug=False):
529 4478 11374.0 2.5 0.1 T_int = self.T.copy()
530 4478 5011.0 1.1 0.0 K = [0.]
531 4478 17374.0 3.9 0.2 pas_de_temps = np.array([0, 0.5, 0.5, 1.])
532 22390 34396.0 1.5 0.3 for h in pas_de_temps:
533 17912 346440.0 19.3 3.4 markers_int = self.markers.copy()
534 17912 179447.0 10.0 1.8 markers_int.shift(self.phy_prop.v * h * self.dt)
535 17912 2151399.0 120.1 21.3 temp_I = markers_int.indicatrice_liquide(self.num_prop.x)
536 17912 72496.0 4.0 0.7 T = T_int + h * self.dt * K[-1]
537 17912 36542.0 2.0 0.4 int_div_T_u = -1 / (self.phy_prop.dS * self.num_prop.dx) * \
538 17912 145813.0 8.1 1.4 integrale_volume_div(T, self.phy_prop.v * np.ones((T.shape[0] + 1,)), I=temp_I,
539 17912 2977345.0 166.2 29.4 dS=self.phy_prop.dS, schema=self.num_prop.schema)
540 17912 127346.0 7.1 1.3 Lda_h = 1. / (temp_I / self.phy_prop.lda1 + (1. - temp_I) / self.phy_prop.lda2)
541 17912 86600.0 4.8 0.9 rho_cp_inv_h = temp_I / self.phy_prop.rho_cp1 + (1. - temp_I) / self.phy_prop.rho_cp2
542 17912 31650.0 1.8 0.3 div_lda_grad_T = 1 / (self.phy_prop.dS * self.num_prop.dx) * \
543 17912 187282.0 10.5 1.9 integrale_volume_div(Lda_h, grad(T, dx=self.num_prop.dx), I=temp_I,
544 17912 2860609.0 159.7 28.3 dS=self.phy_prop.dS, schema=self.num_prop.schema)
545 17912 58233.0 3.3 0.6 rho_cp_inv_int_div_lda_grad_T = self.phy_prop.diff * rho_cp_inv_h * div_lda_grad_T
546 17912 33978.0 1.9 0.3 K.append(int_div_T_u + rho_cp_inv_int_div_lda_grad_T)
547 17912 16489.0 0.9 0.2 if debug:
548 plt.figure('sous-pas de temps %f' % (len(K) - 2))
549 plt.plot(self.num_prop.x_f, interpolate_form_center_to_face_weno(Lda_h) * grad(T, dx=self.num_prop.dx),
550 label='lda_h grad T, time = %f' % self.time)
551 plt.plot(self.num_prop.x, rho_cp_inv_h, label='rho_cp_inv_h, time = %f' % self.time)
552 plt.plot(self.num_prop.x, div_lda_grad_T, label='div_lda_grad_T, time = %f' % self.time)
553 maxi = max(np.max(div_lda_grad_T), np.max(rho_cp_inv_int_div_lda_grad_T), np.max(rho_cp_inv_h))
554 mini = min(np.min(div_lda_grad_T), np.min(rho_cp_inv_int_div_lda_grad_T), np.min(rho_cp_inv_h))
555 for markers in self.markers():
556 plt.plot([markers[0] + self.phy_prop.v * h * self.dt] * 2, [mini, maxi], '--')
557 plt.plot([markers[1] + self.phy_prop.v * h * self.dt] * 2, [mini, maxi], '--')
558 plt.xticks(self.num_prop.x_f)
559 plt.grid(b=True, which='major')
560 plt.legend()
561 4478 17144.0 3.8 0.2 coeff = np.array([1. / 6, 1 / 3., 1 / 3., 1. / 6])
562 4478 116848.0 26.1 1.2 self.T += np.sum(self.dt * coeff * np.array(K[1:]).T, axis=-1)
563 4478 596414.0 133.2 5.9 self.update_markers()
Total time: 27.6902 s
File: <ipython-input-4-6361008d467a>
Function: compare_energy_forme at line 1
Line # Hits Time Per Hit % Time Line Contents
==============================================================
t_fin = 0.2
fig1, ax1 = plt.subplots(1)
for dx in [2*10**-5, 4*10**-5, 7*10**-5, 10*10**-5]:
print('~~~~~~~~~~~~~~~~~~~~~~~~~~')
a = Plotter('decale')
num_prop = NumericalProperties(dx=dx, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy_m
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
l = ax1.plot(t, e, label=prob_ref.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0 # on a mult
print('dE/dt = %g' % dedt_adim)
num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = ProblemConserv2(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy_m
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
l = ax1.plot(t, e, label=prob_ref.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt
print('dE/dt = %g' % dedt_adim)
a.ax.set_ylim(0.,1.05)
a.ax.set_xlim(0., phy_prop.Delta/4.)
ax1.legend()
ax1.grid(b=True)
~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 1.8172585062055175e-06 Cas : convection, rk4, weno, dx = 2.002e-05, cfl = 0.0181544 ========================== dE/dt = -1.69281e-13 dt fourier 1.8172585062055175e-06 Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 2.002e-05, cfl = 0.0181544 ========================== dE/dt = 5.65927e-09 ~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 7.283608525474247e-06 Cas : convection, rk4, weno, dx = 4.00802e-05, cfl = 0.0363452 ========================== dE/dt = -2.01575e-12 dt fourier 7.283608525474247e-06 Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 4.00802e-05, cfl = 0.0363452 ========================== dE/dt = 3.91957e-08 ~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 2.232841866976439e-05 Cas : convection, rk4, weno, dx = 7.01754e-05, cfl = 0.063636 ========================== dE/dt = -8.5763e-12 dt fourier 2.232841866976439e-05 Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 7.01754e-05, cfl = 0.063636 ========================== dE/dt = 1.82845e-07 ~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 4.534064516129032e-05 Cas : convection, rk4, weno, dx = 0.0001, cfl = 0.0906813 ========================== dE/dt = -1.45701e-11 dt fourier 4.534064516129032e-05 Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 0.0001, cfl = 0.0906813 ========================== dE/dt = 4.96835e-07
markers.shift(10**-6/7.)
%matplotlib inline
rc('figure', figsize=(10,7))
rc('figure', dpi=100)
n = 1000
t_fin = 10**-5
fig1, ax1 = plt.subplots(1)
for dx in [5*10**-6, 1*10**-6, 5*10**-7]:
print('~~~~~~~~~~~~~~~~~~~~~~~~~~')
a = Plotter('decale')
num_prop = NumericalProperties(dx=dx, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers=markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(n=n, t_fin=t_fin, number_of_plots=1, plotter=a)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0 # on a mult
print('dE/dt = %g' % dedt_adim)
num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = ProblemConserv2(get_T_creneau, markers=markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, n=n, number_of_plots=1, plotter=a)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt
print('dE/dt = %g' % dedt_adim)
a.ax.set_ylim(0.9,1.05)
a.ax.set_xlim(prob_ref.markers.markers[0][0] - prob_ref.markers.diam / 4., prob_ref.markers.markers[0][0] + prob_ref.markers.diam / 2.,)
ax1.legend()
ax1.grid(b=True, which='both')
~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 1.1335161290322582e-07 Cas : mixte, rk4, weno, dx = 5e-06, dt = 1.13352e-07, cfl = 0.00453406 ========================== dE/dt = -2.04651e-05 dt fourier 1.1335161290322582e-07 Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 5e-06, dt = 1.13352e-07, cfl = 0.00453406 ========================== dE/dt = 2.13251e-07 ~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 4.534064516129032e-09 Cas : mixte, rk4, weno, dx = 1e-06, dt = 4.53406e-09, cfl = 0.000906813 ========================== dE/dt = -3.07688e-08 dt fourier 4.534064516129032e-09 Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 1e-06, dt = 4.53406e-09, cfl = 0.000906813 ========================== dE/dt = 6.01558e-09 ~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 1.133516129032258e-09 Cas : mixte, rk4, weno, dx = 5e-07, dt = 1.13352e-09, cfl = 0.000453406 ========================== dE/dt = -8.76081e-09 dt fourier 1.133516129032258e-09 Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 5e-07, dt = 1.13352e-09, cfl = 0.000453406 ========================== dE/dt = 7.62287e-10